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Abstract. Gravitational- wave signals from inspirals of binary compact objects (black 
holes and neutron stars) are primary targets of the ongoing searches by ground-based 
gravitational- wave (GW) interferometers (LIGO, Virgo, and GEO-600). We present 
parameter-estimation results from our Markov-chain Monte-Carlo code SPINspiral 
on signals from binaries with precessing spins. Two data sets are created by injecting 
simulated GW signals into either synthetic Gaussian noise or into LIGO detector data. 
We compute the 15-dimensional probability-density functions (PDFs) for both data 
sets, as well as for a data set containing LIGO data with a known, loud artefact 
( "glitch" ) . We show that the analysis of the signal in detector noise yields accuracies 
similar to those obtained using simulated Gaussian noise. We also find that while 
the Markov chains from the glitch do not converge, the PDFs would look consistent 
with a GW signal present in the data. While our parameter-estimation results are 
encouraging, further investigations into how to differentiate an actual GW signal from 
noise are necessary. 
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1. Introduction 

Among the sources of gravitational waves (GWs), inspiralling binary systems of 
compact objects, neutron stars (NSs) and/or black holes (BHs) in the mass range 
~ 1 Mq — 100 Mq stand out as likely to be detected and relatively easy to model. For 
ground-based laser interferometers currently in operation (Cutler & Thorne 2002), LIGO 
(Abbott et al. 2009a), Virgo (Acernese et al. 2008) and GEO-600 (Willke et al. 2004), 
the current detection-rate estimates for BH-NS binaries range from 2 x 10~^ to 0.2 yr~^ 
for first-generation instruments {e.g. O'Shaughnessy et al. 2008, Abadie et al. 2010). 
Although the estimates are quite uncertain, detection rates are expected to increase with 
the upgrade to Enhanced LIGO/ Virgo, up to ~ 40yr~^ with Advanced LIGO/ Virgo. 

The detection of a gravitational-wave event is challenging and will be a rewarding 
achievement by itself. After such a detection, measurement of source properties holds 
major promise for improving our astrophysical understanding and requires reliable 
methods for parameter estimation. This is a complicated problem, because of the large 
number of parameters (15 for spinning compact objects in a quasi- circular orbit) and the 
degeneracies between them (Raymond et al. 2009), the significant amount of structure 
in the parameter space, and the particularities of the detector noise. 

In this paper we use an example to illustrate the capabilities of our Markov- 
chain Monte-Carlo (MCMC) algorithm SPINSPIRAL (Van der Sluys et al. 2008a) for 
parameter estimation of binary inspirals with two spinning components, using ground- 
based GW interferometers. In these proceedings we focus on the effects of using 
LIGO detector data versus synthetic Gaussian noise. Earlier studies {e.g. Jaranowski 
& Krolak 1994, Cutler & Flanagan 1994, Poisson & Will 1995, Van den Broeck & 
Sengupta 2007) computed the potential accuracy of parameter estimation {e.g. by using 
the Fisher matrix), but without performing a parameter estimation in practice. Also, 
Rover et al. (2006, 2007), Veitch & Vecchio (20086, 2008a, 2009) explored parameter 
estimation for binaries without spins, described by nine parameters. 

We present the gravitational- wave template used for this study in section [2], and the 
Bayesian framework we employ here in section [31 In section 14.11 we describe the three 
data sets that we analyse in this study; a simulated GW signal injected into synthetic 
Gaussian noise, a GW signal injected into LIGO detector data and a raw LIGO data 
set containing a known artefact of terrestrial origin ("glitch"). We describe the details 
of the MCMC simulations in section 14.21 The analyses of the first two data sets are 
compared in section 14. 3[ and we present our results on the glitch in section 14. 4[ 

2. Gravitational- wave signal and observables 

We analyse the signal produced during the inspiral phase of two compact objects of 
masses M12 in quasi-circular orbit. We focus on a black-hole binary system with 
Ml = 10 Mq and M2 = 1.4 Mq, where unlike in some of our previous studies {e.g. 
Van der Sluys et al. 20086), we do not ignore the second spin to explore the single spin 
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approximation. During the orbital inspiral, the general-relativistic spin-orbit and spin- 
spin couphng (dragging of inertial frames) cause the binary's orbital plane to precess and 
introduce amplitude and phase modulations of the observed gravitational-wave signal 
(Apostolatos et al. 1994). 

A circular binary inspiral with both compact objects spinning is described by a 
15-dimensional parameter vector A G A. Our choice of independent parameters with 
respect to a fixed geocentric coordinate system is: 

A = {A^, rj, log^L, tc, 4>c, tt, cos(5, sin l, ip, 

flspinl , COS 6'spinl , 0spinl , '3spin2 , COS ^spin2 , 0spin2 } , ( 1 ) 

where Ai = and rj = are the chirp mass and symmetric mass 

ratio, respectively; c^l is the luminosity distance to the source; 0c is an integration 
constant that specifies the GW phase at the time of coalescence t^ defined with respect 
to the centre of the Earth; a (right ascension) and 5 (declination) identify the source 
position in the sky; l defines the inclination of the binary with respect to the line of 
sight; and ip is the polarisation angle of the waveform. The spins are specified by 
< aspini.2 = 81^2/ Mf 2 < 1 as the dimensionless spin magnitude, and the angles 
^spini,250spini,2 for their orientations. 

Given a network comprising rzdot detectors, the data collected at the a— th 
instrument (a = 1, . . . , ridot) is given by Xa(t) = Uaif) + ha{t;X), where ha(t; \) = 
Fa^+(t,a,6,ip) ha,+ {t; A) + Fa^y.(t,a,6,ip) ha^y.{t; A) is the GW strain at the detector (see 
Eqs. 2-5 in Apostolatos et al. 1994) and na{t) is the detector noise. The astrophysical 
signal is given by the linear combination of the two independent polarisations ha,+{t; A) 
and ha,x(t; A) weighted by the antenna beam patterns -Fa,+ (^, «, V") and Fa^^it, a, 6, ip). 

The waveform we use includes terms up to 3.5-post-Newtonian (pN) order in phase 
and uses Newtonian amplitudes, with spin effects up to 2.5-pN in phase. We generate the 
waveform templates using the routine LALGeneratelnspiralO with the approximant 
SpinTaylor from the injection package in the LSC Algorithm Library (LAL) (LIGO 
Scientific Collaboration 2007), which closely follows the first section of Buonanno et al. 
(2003). 

3. Parameter estimation: Methods 

In our Bayesian analysis we use MCMC methods to determine the multi- dimensional 
posterior probability-density function (PDF) of the unknown parameter vector A in 
equation [H given the data sets Xa collected by a network of ndet detectors, a model M 
of the waveform and the prior p{X) on the parameters. Our priors are uniform in the 
parameters of Eq.[T](see Van der Sluys et al. (2008a) for details). One can compute the 
probability density via Bayes' theorem 
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where 

C^p{Xa\\,M) OC exp Xa\ha{\) > "3 < ha{\)\ha{\) >^ (3) 

is the likelihood function, which measures how well the data fits the model M for the 
parameter vector A. The term p{xa\M) is the marginal likelihood or evidence. In the 
previous equation 

is the overlap of signals x and y, x{f) is the Fourier transform of x{t), and Sa{f) is the 
noise power-spectral density in detector a. The likelihood computed for the injection 
parameters Cinj = p{xa\\mi, M) is then a random variable that depends on the particular 
noise realisation Ua in the data Xa = ^(Ainj) + na- The injection parameters are the 
parameters of the waveform template added to the noise. We define the signal-to-noise 
ratio (SNR) of the injection to be: 

SNR = . < ^l^(^-^)_^> . (5) 

^J< /l(Ainj)|/l(Amj) > 

From here on, we use the expected value of the SNR, which is equal to the square root 
of twice the expectation value of log£inj: 



SNR=V</^(Ai„j)|MAinj)>. (6) 
To combine observations from a network of detectors with uncorrelated noise 
realisations (this is the case in this paper as we use two non-co-located detectors) we 
have the likelihood p{x\X, M) = Y[a=i p{^a\X, M) , for x = {xa : a = 1, . . . , n^ct} and 

The numerical computation of the PDF involves the evaluation of a large multi- 
modal, multi-dimensional integral. Markov-chain Monte-Carlo (MCMC) methods {e.g. 
Gilks et al. 1996, Gelman et al. 1997, and references therein) have proved to be especially 
effective in tackling this numerical problem. We developed an adaptive (see Figueiredo & 
Jain 2002, Atchade & Rosenthal 2005) MCMC algorithm to explore the parameter space 
A efficiently while requiring the least amount of tuning for the specific signal analysed; 
the code is an extension of the one developed by some of the authors to explore MCMC 
methods for binaries without spin (Rover et al. 2006, Rover et al. 2007). We implemented 
parallel tempering (Hukushima & Nemoto 1996, Hansmann 1997, Rover 2007) to 
improve the sampling. It consists of running several MCMC chains in parallel, each 
with a different "temperature", which can swap parameters under certain conditions. 
Only the T = 1 chain is currently used for post-processing. 

In Eq.[7]we applied Bayes' theorem to obtain the probability of a specific parameter 
vector value (A) given the observed data x and the model M. The theorem can also be 
applied to compute the probability of a specific model Mi given the observed data: 

p(Mi)p(x\Mi) 
p{x) 
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We compare the two models Mj and Mj by computing the odds ratio: 

^ p{Mi\x) ^ p{Mi)p{x\M,) ^ pjMj) 
"'^ p{Mj\x) p{Mj)p{x\Mj) p{Mj) ^ ' 

where 

p{xm ...1 

is the Bayes factor of the two models, and we recognise the evidence p{x\Mi) from 
Eq.[71 The evidence must be marginalised over the parameters of the model in order to 
compute the Bayes factor: 

p{x\Mi)= I p(A|M,)p(x|A,Mi)dA. (11) 

There are existing algorithms dedicated to the computation of this integral, and of the 
Bayes factor. For instance, nested sampling (Skilling 2006) has been shown to be very 
efficient in the case of non-spinning gravitational- wave sources (Veitch & Vecchio 2009), 
and can in addition be used to produce PDFs of the parameters. As a by-product of the 
exploration of the parameter space with MCMC, it is possible to compute the evidences 
of the models used. We have implemented the harmonic-mean method (Newton & 
Raftery 1994), in which the evidence is approximated by: 

N 

p{x\Mi) ^ Y.PiMM^)p{x\\k,Mi)V-^^, (12) 

k=l 

where {A^ : = 1, . . . , A^} is the set of points sampled by the MCMC, and V^^ is the 
volume of parameter space associated with the point A^. Since the MCMC algorithm 
samples according to the posterior (and, up to a proportionality constant, converges 
towards posterior PDF), the density of points in the chain at a certain location A^ in 
the parameter space A will become proportional to the posterior for large A^. It follows 
that 

lim = , (13) 

^'^ p(Afc|M,)p(f|Afc,M,) 

with Ui a proportionality constant. We then have p{x\Mi) ^ J2k=i'^i = N ai, and 
obtain the estimate for ai by considering the whole parameter space volume Vt: 

Vt^Tv,=T^ . (14) 

ti ' tip{Xkmp{x\Xk,M,) 



Finally, 



p{x\Mi) ^ NVt 



N 

E 



1 



(15) 



[t^ip{Xk\Mi)p{x\Xk,M,) 

which is the harmonic mean of the posterior values sampled by the MCMC. The issue 
with this method is that it gives too much weight to low-posterior points, which lie 
in a part of the parameter space that is badly sampled, by design, by the MCMC. 
The estimate of the evidence is then very sensitive to the quality of the sampling of a 
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particular run. We are looking into other algorithms in order to remedy this problem, 
e.g. by using the higher-temperature chains produced by parallel tempering (Earl & 
Deem 2005) (we currently use the T = 1 chain only), or by using a well sampled subset 
of points (Van Haasteren 2009) to estimate the probability constant Oj. A summary of 
the methods used in our MCMC code was published in Van der Sluys et al. (2008a); a 
more complete technical description of the SPINSPIRAL code will be available in (Van 
der Sluys & al. in preparation). 

4. Parameter estimation: Results 

4.1. Data sets 

For these proceedings, we analyse three different data sets, each containing the data for 
the 4-km LIGO detectors at Hanford (HI) and Livingston (LI): 

DSl: a coherent software injection with a total SNR of 11.3 into synthetic Gaussian, 
stationary noise, simulated for the HI and LI detectors; 

DS2: a coherent software injection of the same signal, with a total SNR of 11.3, into 
"quiet" LIGO detector data from HI and LI; 

DS3: raw LIGO data from HI and LI, containing a known, coincident glitch of seismic 
origin, with a total SNR of 11.3. 

For the data sets DSl and DS2, the injected signal is that of a 10 spinning BH 
and a 1.4 Mq spinning NS in an inspiralling binary system. A low-mass Compact Binary 
Coalescence Group search (Abbott et al. 20095) does not produce a GW trigger for the 
data segment DS2; hence we designate it "quiet". The distance of each of the injections 
is scaled to obtain an SNR of 11.3, equal to that of the glitch in DS3, but computed 
with different waveforms: a SpinTaylor waveform (see section[2]) for DSl and DS2, and 
a non-spinning, 2-pN waveform (see section l^^ for DS3. The other parameters of the 
injection are: 

\ = {M = 2.99 Mq, ri = 0.107, du tc, 0c = 85.9°, a = 17.4 h,6 = 61.6°, 
^ = 52.8°, = 11.6°,a,pi,i = 0.6,^spini = 78.5°,0spini = 63.0°, 

aspin2 = 0.4, ^spin2 = 120.0°, 0spin2 = 315.1°}, (16) 

where we assigned a spin of 0.4 to the neutron star, which is higher than astrophysically 
plausible, for testing purposes only. In DS3, no signal is injected. For our analyses, we 
use the data of both 4-km LIGO detectors HI and LI. 

4.2. MCMC simulations 

The MCMC analysis that we carry out on each data set consists of 10 independent 
Markov chains, each with a length of about a million iterations and composed of 5 
chains at different temperatures for parallel tempering. From now on, we will refer 
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to the T = 1 chain as the chain, since the hotter chains were not used in the post- 
processing. The part of the chains that is analysed is that after the burn-in period (see 
e.g. Gilks et ah 1996), the length of which is determined automatically as follows: we 
determine the absolute maximum likelihood log(£niax)5 defined as the highest value for 
log[p(a;|A, M)] obtained over the ensemble of parameter sets A in any of our individual 
Markov chains. Then for each chain we include all the iterations after the chain reaches 
a likelihood value of log(£max) — 2 for the first time. This results in a convergence 
test as well, since some of the independent chains may not reach this threshold value. 
Typically, we demand that more than 50% of our chains meet this condition before 
we consider the MCMC run as converged, although we consider results as robust if 
they have a convergence rate of 80% or more. This convergence test is a measure of 
the quality of our sampling in a given number of iterations. All our Markov chains 
start at values that are randomly offset from the injection values. The starting values 
for Ai and tc are drawn from a Gaussian distribution centred on the injection value, 
with a standard deviation of 0.025 Mq and 10 ms respectively. In real analysis, the two 
Gaussian distributions are centred on the values from the template bank based search of 
the Compact Binary Coalescence group (Abbott et al. 20096) which will have triggered 
the MCMC followup. The other thirteen parameters are drawn uniformly from their 
allowed ranges. SPINSPIRAL needs to run for typically a few days in order to show the 
first results and a week or two to accumulate a sufficient number of iterations for good 
statistics, each chain using a single 2.8 GHz CPU. 

4.3. Analysis of data sets DSl and DS2 

We analysed the data sets DSl and DS2 as described in section 14.11 and the results 
of both analyses passed the convergence test described in section 14.21 with convergence 
rates of 70% and 80%, respectively. The resulting one-dimensional marginalised PDFs 
from both analyses are shown in figure [H 

Table [H shows the median and the width of the 95%-probability ranges for each 
parameter. The differences we find between the results for DSl and DS2 may be 
attributed to the particular noise realisations in this example, and most parameters 
yield similar PDFs and accuracies. 

The PDFs of the parameters that describe the spin of the NS follow the prior 
distributions in both runs. This justifies ignoring the NS spin (by fixing aspin2 to 0.0 in 
the recovery template) for this mass ratio (Van der Sluys et al. 20086). For each of the 
two data sets, DSl and DS2, we computed the Bayes factor to compare the evidence 
for the following two models: Mi. a 3.5-pN inspiral waveform embedded in Gaussian 
noise, and M2: Gaussian noise only. The values are listed in table |2j In both cases, 
the Bayes factor is large, providing strong evidence for a GW signal in the data. The 
difference in Bayes factor between DSl and DS2 is attributed to an inherent spread 
due to different noise realisations, and the uncertainties of our method to estimate the 
Bayes factor (section [3]). The results in this section show an illustrative example, but 
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Figure 1: One- dimensional marginalised PDFs for all 15 parameters from our analysis 
of data sets DSl (hatched upward; red in the online colour version) and DS2 (hatched 
downward; blue in the online colour version). The vertical dashed lines mark the 
injection values. 

cannot be used to draw firm conclusions. However, it is clear that they warrant a larger, 
systematic study of these phenomena with the methods described here. 

4.4- Analysis of data sets DS2 and DS3 

On November 2nd 2006, seismic activity at Hanford and Livingston resulted in a 
coincident "glitch" in the data from the HI and LI LIGO detectors. These glitches 
were recovered by the Compact Binary Coalescence detection pipeline at an SNR of 11.3, 
using non-spinning, stationary-phase-approximation templates, Newtonian in amplitude 
and 2.0-pN in phase (Abbott et al. 20095). We defined the corresponding data set as DS3 
in section 14.11 and analysed the data as if it had yielded a GW trigger. The convergence 
test from section 14.21 yields a 20% convergence rate, which results in our rejection of 
the results as not converged. However, when we nevertheless construct the marginalised 
one-dimensional PDFs from the data of the two converged chains (because of the small 
number of data points, the resulting PDFs may not be very accurate), they are similar 
in appearance to those from DS2 (see figure [2]). The Bayes factors in table |2] even 
suggest that the data set DS3 is more consistent with containing a GW signal than DS2 



Parameter estimation of gravitational-wave signals with spin effects 



9 



Table 1: Median and width of the 95%-probabihty ranges for each parameter of the 
analyses of data sets DSl and DS2. The column recovered indicates whether or not the 
95% range includes the injection value. 







DSl (synthetic 


noise) 


DS2 (detector 


noise) 




injection 


median 


95% width 


recovered 


median 


95% width 


recovered 
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3.006 


0.294 
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3.041 
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yes 
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yes 
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21.240 
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24.144 


17.238 


yes 
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-0.013 


0.024 
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0.006 


0.019 


yes 
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342 398 


ye& 
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343.175 




a{h) 


17.380 


11.684 


5.349 


no 


17.786 


6.320 


yes 


sn 


61.642 


49.326 


64.346 


yes 


58.390 


39.796 


yes 


HI 


52.753 


67.056 


110.735 


yes 


46.850 


122.787 


yes 




11.459 


93.162 


176.358 


yes 


88.706 


173.869 


yes 


^spinl 


0.600 


0.658 


0.594 


yes 


0.804 


0.478 


yes 


^'spinl (°) 


78.463 


85.490 


83.110 


yes 


89.225 


85.787 


yes 


^spinl (°) 


63.025 


57.171 


335.592 


yes 


263.014 


345.700 


yes 


^spin2 


0.400 


0.532 


0.945 


yes 


0.475 


0.940 


yes 


^'spin2 (°) 


120.000 


94.687 


150.544 


yes 


89.406 


146.101 


yes 


0spin2 (°) 


315.127 


181.959 


327.603 


yes 


184.681 


339.071 


yes 


Afi (Mo) 


10.002 


8.533 


8.849 


yes 


6.421 


6.536 


yes 


Ah (Mq) 


1.400 


1.598 


1.277 


yes 


2.036 


1.564 


yes 



Table 2: Bayes factors B12 between the models Mi: a 3.5-pN inspiral waveform 
embedded in Gaussian noise, and M2: Gaussian noise only (section l4.3p for data sets 
DSl and DS2 (see section E]). 



DSl (Gaussian noise) DS2 (detector data) DS3 (glitch) 
loggBi,2 52.9 43.5 68.5 



(with the caveat that the SNRs of DS2 and DS3 were not computed the same way). On 
the other hand, the low value for the median of rj (0.05) corresponds to a mass ratio of 
18, which is near the limit of the regime where post-Newtonian expansions are valid. In 
particular, a small value for eta suggests a slow frequency evolution which may indicate 
a spike in the frequency spectrum that dominates the signal. In addition, we find that 
the sky map for DS3 does not display the (parts of a) sky ring that is expected for 
an analysis using two non-co-located detectors (see e.g. Raymond et al. 2009). These 
results indicate that we should thoroughly verify our tests, such as the convergence 
criterion described here, using a large number of different glitches. 
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Figure 2: One-dimensional marginalised PDFs of a few selected parameters from our 
analysis of data set DS3. The vertical dashed lines indicate the median of each PDF. 

5. Conclusions 

We have developed the code SPINSPIRAL which can do a complete parameter analysis 
of the gravitational-wave signals from quasi-circular compact-binary inspirals. We 
presented an example of the analysis of software injections into both simulated Gaussian 
noise (DSl) and LIGO-detector data (DS2). We also presented an analysis of a data set 
containing no injection, but a "glitch" coincident in two LIGO interferometers (DS3). 
These examples demonstrate a remarkable similarity between the results obtained from 
a GW signal injected in Gaussian noise and a similar signal in detector data. The Bayes 
factors are also similar, where we note that our present technique for computing the 
Bayes factor yields estimates with significant variance, and more precise estimates should 
be possible in the future. In addition, we find that although the Markov chains in the 
analysis of a coincident glitch in LIGO data do not converge, the resulting PDFs could 
look remarkably consistent with a simulated GW signal. We plan to run our code on a 
very large number of coincident triggers from the LIGO Compact Binary Coalescence 
search pipeline (noise events that are somehow being registered as resembling a binary 
inspiral) in order to get a good sense of how to distinguish them from actual inspirals. 
We conclude that further, detailed investigations are necessary to ensure we can rely on 
the robustness of our tests. 
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